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INTRODUCTION 


There has recently been an increased awareness of the fact that pollutants which 
are normally invisible to the naked eye may pose health hazards. Thus, major efforts 
have been mounted to control various forms of those hazardous pollutants. 

For effective control, effective prediction is necessary. The forecasting of air 
pollution concentrations presents a very difficult problem of predictibility since many 
variables influence the concentration of pollutants that are measured at receptor sites. 
Numerous examples of such variables can be given, such as source strength, source geo- 
metry, meteorological influences, chemical reactions involving the pollutants, etc. 

The case of automotive pollution emanating from a line source is even more 
complex than the case of a point source such as a factory. It has been remarked in a 
recent NASA publication [1] that air pollution models which deal with diffusion from 
highways and airports are not as refined as those that deal with general urban models. 
The purpose of this document is to produce a more reasonable model dealing with 
highways. 

The method that will be followed in the derivation of applicable equations is 
analytic rather than purely numerical. A sophisticated and satisfying treatment could be 
made via numerical integration of the diffusion equation, but the computational time 
involved in such a scheme would be quite high. Additionally, numerical integrations on 
systems of this sort tend toward instability. 

The analytical model which is derived here is a compromise between the simple 
approximations such as occur in Reference 2 and numerical methods. 

The following flow chart (Fig. 1) illustrates the structure of the overall model. 



Figure 1. General Flow Chart. 









An expanded structure for the block labeled “diffusion model” is required, and we 
consider inputs to such a model. The basic diffusion model which governs the flow of 
pollutants is fundamental to the entire discussion. The mechanics of the nondiffusing 
air is accounted for via meteorological inputs. The source term can be considered to be 
composed of two components, namely road geometry and traffic density. The important 
area of chemical reactions of the pollutants will not be considered at this time. The 
following flow chart (Fig. 2) summarizes the details of the internal structure for the 
blocks labeled “meteorological conditions” and “roadside pollution concentration.” 



MODEL 


Figure 2. Boundary conditions. 

The above outline prescribes an approximate diffusion model which is coupled to 
both the prevailing meteorological conditions and the (time variable) traffic density. The 
results of the study provide a rapid method of estimating 'the concentrations of pollutants 
from a line source. Based upon the calculations given here, a full scale simulation 
involving direct numerical integration of the applicable partial differential equations could 
be undertaken. 

In the development, a general separation of the diffusion equation will be accom- 
plished via a product assumption. This requires several subsidiary physical hypotheses, 
such as thermodynamic independence of the diffusion tensor and averaging of the winds 
in space and time. The separates chosen are of such form that boundary condition 
matching is facilitated. The specified separation is followed by boundary condition 
matching on meteorological parameters and roadside pollution concentration. Finally, 
the varying road geometry is taken into account by eithwi numerical integration or 
segmentation or an approximate analytic method. The theory is applied to yield parametric 
results and the numerical behavior of the equations is demonstrated. 
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THE DIFFUSION EQUATION 


The standard equation for a diffusing gas [2] is 

+ ^ • Vx = V- (K Vx) , (1) • 

dt 

where 

X is the concentration of the diffusing gas, a function of both time and space variables, 

V is the vectorial wind velocity, which will be taken as an average value, 

K is the diffusion tensor, 

t is the time, and 

V is the spacial gradient operator. 


It is possible, via a rotation to diagonalize K so that 

AA AA A A 

K = kjj i i + ky j j + k k , (2) 

where t, and Ic represent cartesian unit vectors. The components of K (kj^ , ky , 
and k^) represent the so-called eddy diffusivities, and are functions of the down-wind 

distance. The physical basis of a variable diffusivity is that the pollution cloud spreads 
as it leaves the road and interacts with eddies which are of the scale of the cloud. Thus, 
taking k^ , ky , and k^ as constants is not physically satisfactory. By Reference 2, we 
take 

k^ = 0 (by superposition), 

ky = ky(x) , and (3) 

k^ = k^fx) 
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u. 



Assuming that the vector velocity, V , lies along the x axis yields 


—A 

V = u i 


where u = IVl . (The case where V does not lie along the x axis will be considered 
at a later time.) 

Inserting (2), (3), and (4) into (1) yields 


ax 3x a'x _ , , 9"x 


This equation is at least partially separable. Assume that 


X = T(t) Y(y) W(x,z) 


so that (5) becomes 


r ^ ^ 

T w 3x 


Y" k^Cx) g2'V 

kyCx)— + -~ 

y Y w 3z^ 


where the prime (') denotes differentiation with respect to the argument. Isolation of 
the only component dependent upon T yields 


r , Y" 3^W u 3W 

— - - = k.,(x) — + — — - — — 

T y Y w 3z^ w 3x 


where is a separation constant of unknown properties. Thus 


T = A* e 
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where A* is an integration constant. Isolating the Y component of (8) yields 


Y” _ 1 I" 2 ^ 3^1 

Y ky(x) L w 3z^ ^ w 9xJ 

where - is a second separation constant. Then 

Y = B* cos(/3y + 0y) 

and B* , 0y are integration constants. 

The residue of (10) may now be written as 


( 10 ) 


( 11 ) 


l<2(x) 


d^w 




( 12 ) 


It would be physically reasonable to require that u = u(z) and proceed with the inte- 
gration of (12) on that basis. A difficulty arises in that no method of integration has 
been found for this equation with variable and ky unless u is assumed to be con- 
stant, This is, mathematically, why most authors always use average values of Ti. To 
obtain an initial separation it was necessary to assume that u was averaged in time and 
along y. It is now neces.sary to assume that Ti is averaged in z as well. In Reference 
3 it is stated that one almost never has enougli data to establish the dependency of u 
on z in any case. Thus, we assume that 


W(x,z) = X(x) Z(z) 


(13) 


so that (12) becomes 


Z” 

z" 


_J 

kj,(x) 



+ u 


X' 


+ /J' ky(X) 




-7‘ 


(14) 


The total differential equation for Z is, then. 
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/!' + 7 ^ 


0 


( 15 ) 


The solution of equation (15) could be written exactly as (11) was written as a 
solution for (10). But it is convenient to choose a special form with an eye toward 
matching boundary conditions. Thus we write, 


Z 


c* 


oo 

cos 7 (z - H) + cos 7 (z + H) + ^ [cos 7 (z - H - 2nl) 

n=l 


+ cos 7 (z - H + 2nl) + cos 7 (z + H - 2nl) + cos 7(2 + H + 2nl)J j 

(16) 

(The parameters H and 1, arbitrarily introduced into (16), have physical significance. 

It will later be shown that 1 is the depth of a mixing layer, i.e., an upper lid which 
is not penetrated by the pollution cloud. H is the effective height of emission of 
the gases which arise due to the fact that such gases are normally warmer than the 
ambient air.) 

We also have, from (14), 


X' = 4 [- 7 ' k^(x) + ky(x)J X 


(17) 


which has the solution 


X = D* exp 


X 


?2 . X 


- — / k^(x) dx f ky(x) dx + 4 - (x - x^) 

u ^ u “'v ^ u “ 


(18) 


where Xq is an arbitrary reference length. (Later Xq will be chosen to be the half- 
width of the road.) 

For algebraic convenience, we define 
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X 


, and 


f(x) =J k^Cx) dx 


( 19 ) 


h(x) = J ky(x) 


dx 


Using (19) in (18) and the resultant equation, along with (16), (13), (11), and (9) in (6) 
gives 


X = A* B* C* D* e 


V VI 7 u 

e V ' e. 


f(x) 0^h(x) 


u 


cosOJy + ^y) < cos y(z - H) 


+ cos yiz + H) + ^ jcos y(z - H - 2nl) + cos y{z - H + 2nl) 
n=l 


+ cos y(z + H - 2nl) + cos y(z + 
which is the formal separated solution of (6). 


H + 2nl)J I 


( 20 ) 


BOUNDARY CONDITION MATCHING 


Since the separated equation is degenerate in the (t,x) pair, it is obvious that 
boundary conditions can be matched only on z , y , and (t,x). We begin with z . The 
eigenvalues, y , form a continuous spectrum so we can write 


7^f(x) 


= xtyJ 


COS 


oo 

7(z - H) + cos 7(z + H) + ^ jcos 7(z - H 

n=l '■ 


- 2nl) 


+ cos y(z - H + 2nl) + cos 7(z + H + 2nl) + cos 7(z + H - 2nl)J 


d7 (21) 


From Reference 5, eq. 3.896 we have 
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( 22 ) 



/ 


e ^ cosq2(x + X)dx = e cos qj 

Qi 


X 


Identification of variables between (22) and (21) yields 


oo 7^f(x) 

f ^ ^ cos 7(z ± H ± 2nl) d7 

— OO 



(z ± H ± 2nl)2 u 
4f(x) 


(23) 


Similarly, (21) can be written as 


X 


OO 

= XTZf 


— OO 


/?"h(x) . 

e “ cos(/3y + 0y) dy 




cos 0y 


(24) 


Thus, including all possible values of j3 and 7 shows that (20) may be written, with 
A = A* B* C* D* cos 0y , as 


Aim 



-^(z-H)2 _^(z + H)2 

e + e 


00 

E 

11= 1 


-■^(z-H-2nI)2 


+ e 


- ^ (z + H + 2nl)2 


+ e 


- ^ (z + H - 2nl)2 


+ e 


- ^ (z - H + 2nl)2 


(25) 


It is now convenient to modify out notation slightly to be in agreement with the 
standard literature. The quantities u/f and u/h are not generally tabulated. Instead 
we find, as in [3] , direct tabulations of the quantities Oy and Oj,. Physically Oy and 

are quantities such that ~ 97 percent of the mass of the pollutants are contained in 
a cloud of' radius 2.15a. Mathematically, 
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( 26 ) 


^2f 

= + / — 
u 


are given. Inserting Oy and in (25) gives the more compact form 


X = 


2A7T 

Oy O^ 




M 1 (Z - H)^ 1 (Z + H)^ 

2 Oj,^ 2 

e 2 + e ^ 


- S 


n=l 


(z - H - 2nl)2 (z - H + 2nl)2 (z + H + 2nl)^ 


2a,^ 


+ e 


2a^^ 


+ e 


2a^^ 


+ e 


(z + H - 2nl)2 


(27) 


It is now possible to illustrate an important set of properties of the “Z” portion 
of (27). We have 


9Z 

9z 


1 


a 


z 


2 


(z - h) e 


1 (z - H)^ 

2 


+ 


(z + H) e 


1 (z + H)2 

2 


OO 

+ z 


n=l 


(z - H - 2nl)2 
2a 

(z - H - 2nl) e ^ 


(z - H + 2nl)^ 


2a~^ 

+ (z - H + 2nl) e ^ 


(z + H + 2nl)2 


(z + H - 2ni)2 


2o^^ 

+ (z + H + 2nl) e ^ 


2a~^ 

+ (z + H - 2nl) e ^ 


(28) 
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Notice tliat 


3Z 


3z 


z=0 


0 


(29) 


Similarly, an expansion of the series shows that 


3Z 


3z 


z=l 


0 


(30) 


From [4] , equations (29) and (30) ensure that the ground and a mixing layer of 
arbitrary depth, 1, are impervious to pollution transfer. 


Also, if we allow 1 to become large in (27) we find that Z can be written as 
1 (z - H)2 1 (2 + H)2 


Z = XTY 


+ e 


which corresponds to the physically unrealistic case of a superadiabatic atihosphere of 
infinite depth. 

The next modification of (27) to meet the specific model we wish to consider 
involves the fact that (27) represents a point source, not a distributed source such as a 
road. The following figure is applicable. 


X (WIND. DIRECTION) 



(SOURCE POINT) 


Figure 3. Distributed Source Geometry. 
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(-y^, Vq are, respectively, the Eastern and Western terminus point of the road segment 
under consideration.) 


We have, from (27), 


X 


XTZ e 



for a point source at (0,0). But for a point source at (0,y*) of strength Qdy* per unit 
distance, we must have 


(y - y*) ^ 

20y2 


dx = QXTZ e 


dy" 


or 


X = QXTZ / e 
-Vo 


(y - y*)^ 


dy"* 


Setting 


0 


(y* - y) 

\/2ay 


we have 


d© 


dy* 

v/2oy 


And when 


y* = ±Vo 


/ 
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wc must have 


0 = ± 


(Mo + y) 

\/2 Oy 


Thus, 


vo-y 


X = QXTZ 


r yAoy 

— 0* 


1 

0 

+ 

1 

I (yo + y) 

e d0 = y/W Gy QXTZ ( 

erf 

v/TOy 




+ erf 


yp-y 

\/2 Oy 


or 


2 A(7T)^/2 q 

X = e 




erf 


(yo + y> 




+ erf 


yo-y 


y/2 a. 


1 (z- 


{e 


2 a 


1 (z + H)2 

^ oo 

+ e ^ + ', 


(z - H - 2nl)2 (z - H + 2nl-r (z + 


2a,2 


+ e 


2a, 


+ e 


+ e 


(z + H + 2nl)^ 

2^? 


(31) 


H)2 


H - 2nl)2 

2^? 


(32) 
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Equation (32) is not yet normalized and does not yet account for traffic density. 
There would be a certain elegance involved in a car-by-car trace which gives rise to a 
complicated coupling between y and t , but this is hardly necessary. We will, instead, 
simply assume that the travel time of any car on the road is short compared to the time 
over which the pollution cloud is observed. Equivalently, we shall assume that the 
variation of traffic at the center point of the road represents the variation of traffic along 
the entire road, and match our boundary condition in time at this point. 

In order to match the traffic density boundary condition, we shall make use of the 
the remaining separation parameter (a) and the “constant” A. Note that no violation of 
equation (1) results if we assume that A is a function of a. Furthermore, in order to 
obtain a normalized function we evaluate x as the point (x^, 0, 0, t) as 


X(xQ,0,0,t) = 


8 A(o()(7r)^/^ Q -a^t / ^o 

e / erf 


^z(^o) \ s/Tay(Xo) 

(H + 2nl)2 (H - 2nl)2 


1 — 

2 


OO 

* E 


2a, 


+ e 


2a, 


n=l 


For algebraic convenience, we label 


(33) 


8(7t)'*/^Q 
a^/Xo) 


«!> = 


erf 




1 


2 a, 


^ - z 

n=l 


(H + 2nl)2 (H - 2nl)2 


2a, 


+ e 


2a, 


(34) 


so that (33) becomes 


X(xQ,0,0,t) = A(a) e 


-a^t 


(35) 


Accounting for all possible values of a gives 


13 



( 36 ) 


w - 

r -0^ t 

X(XQ,0,0,t) = 4> J A(a) e da 

where the lower limit is deliberately unspecified. 

We wish now to assert that at (x^,0,0,t), the pollution density at roadside is given 
by 

X(xo,0,0,t) = QN(t) , (37) 

where Q is the emission per car per unit distance and N(t) is the number of cars per 
unit time. Thus given N(t) we wish to solve for A(a) such that 


OO 

QN(t) = $/ 


-a^t ^ 

A(a) e da 


(38) 


where N(t) may be expanded into an arbitrary series which is subject only to well 
behaved numerical behavior properties. (The possibility that N(t) could be represented 
by a single, simple, analytic function is so remote as to be negligible.) 

Three cases will be considered. In the first case we shall deal directly with a^ 
and the range of integration will be (-«»,'»). Secondly, if we assume, that a can be 
written in place of a^ , then A(a) becomes the inverse Laplacian of Q N(t)/4> . The 
integration range, in that case, is (0,°®). Finally, if we assume a to be a purely imaginary 
number, which implies that A(a) is the inverse Fourier transform of Q N(t)/^ . 

For the first possibility, assume that 


A(a) = ^ ap(a^)^ 

p=l 


(39) 


so that (38) yields 


QN(t) 



^ 2p -a^ t 

L 3pa e da 

p=l 


^ ^ (2p- 1) W ^ 

*1, “p 


(40) 
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(Reference 5, eq. 3.461). Since N(t) is an empirical function, assume that it may be 
presented as 


M b. 


N(t) = y — ^ 


Then we must have 


2 ” b. 


"P 


(2p - 1 ) ! ! y/ir 4» 


so that 


M 2 P b„ «2P 
A(c) = I ” 


(2p- 1) !!v^ <P 


(41) 


(42) 


(43) 


Finally, the expression for equation (32), can be written as 


X 



M 2 P b„ «2p 
V ___P 

(2p- 1) !!v^ 




(44) 


We next assume that a may have only positive values and choose 


M 

A(oi) = Yj 3p cos p a , (45) 

P=1 
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so that 


^ ^ -at ' ^ ^ 

X(x_, 0 , 0 ,t) = Q N(t) = $ J Yi a cos p a e da = ^ .2 ~„2 

o p=l - p=l " 


( 46 ) 


Assuming that x is readily fit as 


M 

= qZ p 

p=l 


bp t 

+ P' 


we identify 


Qbp 

“p $ 


SO that 


( 47 ) 


( 48 ) 


X 


2 (tt)^/^ Q 





cos p a e 




( 49 ) 


Finally, we choose a purely imaginary separation constant and define 


A(a) = 


M 


Z 

p=i 






( 50 ) 
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riU'U 


M . 


X(xo,0,0,t) = QN(t) = <pf 


V2PZ 


/ 7T \ 

sin ( — — + — I 

Y pp V 


-iat 


M 


p=l 


( 51 ) 


Thus, if we can write, conveniently. 


M 

N(t) = Yj ^^P ^ 

p=l 

we must have 



The choice of equation (44), (49), or (54) is arbitrary from the theoretical point 
of view. Numerically, however, the behavior of the three equations differs markedly. 
Whichever of these best fits the traffic pattern should be used. 
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Kquation (54) may be fully written as 


X(x,y,z,t) = 


Q 

4 a7(x)‘ 


Vv^oy/ X^-Oy/. 


erf 


M 

t <>p“"Pp 

p=l 




< 7 y(Xo) 


(z - H)* (z + H)2 




V 

Li 

n=l 


(z - H - 2nl)2 (z - H + 2nl)* (z + H + 2nl)’ (z + H - 2nl)’ 


2 oz 


2o,2 2o = 

+ e ^ + c * 


+ e 


2a,‘ 


( 55 ) 


H» 


2 »z(x„)’ 


( (H * 2nl)» ( H - 2nl)» \ 
n=l 


Equations (44), (49), and (54) represent three possible solutions out of an infinite 
number of choices which could be made for A(a), and they are convenient from certain 
points of view. An additional purely numerical scheme is also worth considering. These 
equations all share the property that they are degenerate in the variables t and x. This 
is a direct result of the assumption that kj^(x) vanishes. Set, for convenience 


X-Xq 

w(t,x) = t r — ^ 0 

u 


and consider a coordinate X| > x^, and a time t = tj . Then 
3t* ; (t*<t) Aj^w(t*, Xq) = w(tj,X|)J 


where we can explicitly define 



n 


This represents a retarded time wave. The practical result of this observation is that if 
we define the traffic density as a dense function (perhaps even via interpolation), then 
the representation of the traffic via calculation of polynomial coefficients becomes 
unnecessary. 
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ROAD GEOMETRY 


It was earlier stated that we would temporarily assume that the road was per- 
fectly straight and that the vectorial wind was orthogonal to the road. In practice, of 
course, roads are rarely straight and this fact must be included in the theory which has 
been presented. 

Several methods of approach suggest themselves. We could, for example, curve 
fit a polynomial between the x and y coordinates of the road, obtaining a relationship 
such as 


y* = F(x*) 


(56) 


and integrate an equation similar to (31), after replacing x by x - x* and y by 
y - y* . Such a process is possible, but would most likely prove unsatisfactory due to 
the complexities of curve fits within curve fits and the machine time required for com- 
plicated numerical integrals. 

A slightly modified approach would be to break the road into linear segments 
and not assume that the wind is orthogonal to the road. That is, assume that (56) is a 
linear relationship via 


y* = X* tan \j/ 


(57) 


where i// is the angle Ijetween the y axis and the road (Fig. 4). 


X (WIND DIRECTION) 



Figure 4. Nonorthogonal Wind Geometry. 
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riuis, (55) would he wrilteu as 


X(x,y,z,t) = - a^(\^)ay(\Q) J 

-Vo tan \p 


(y - X* tan 4/)^ 

2ay^(x-x*) 
e ” 


M 

Z. '^p Pp 

p=l 


X - X* - X, 


t - 


ay(x - X*) a^(x - X*) 


(z - H)2 (z + H)2 


2a2(x-x*) 2a,^(x-x*) 

^ + e 


OO 

- s 

n=l 


(z - H - 2nl)2 (z - H + 2nl)2 (z + H + 2nl)2 


2a^^(x-x*) 2a^*(x - X*) 


+ e 


+ e 


20^2 (X - X*) 


+ e 


(z + H - 2nl)^ 

lOy^iX - X*) 


dx* 


which is a none too pleasant integral.* 


A third possibility, an approximation that offers some hope of numerical applica- 
tion is as follows. Consider a road of arbitrary geometry. Segment the road into a large 
number of individual straight segments and apply, say, (55) to each of the pieces (Fig. 5) 
Subsequently, the segments are added linearly to obtain the total air pollution. Graph- 
ically we have: 


X 



Figure 5. Nonlinear Source Geometry. 

* An analytical approximation to this integral is developed in Appendix B. 
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To rewrite (55) in such a way as to account for such segmentation, assume that 
the road is broken into k segments, each of length 2L (Fig. 6). Label a generic seg- 
ment by the subscript i , and assume that it is centered at the point xj , yj . 


X 



If we assume an imaginary separation constant, then the new form is available directly 
from (55). Thus we have an accounting for all segments. 




l2 a,(x- )C|) 


(z - H)= 


(y - Yj - L,) (y - Yj + Lj) 

erf + erf 


\/^ tJy(X — Xj) Qy(y>. Xj) 

\/2 0yU^)) - \N/2Cy(X„y 


Yj •’p ’’p 

p=l 




(z + ny 


20 y/(x-Xj) - Xj) 

e + e + ^ 

n=l 


(z + H + 2n\y (z + H - 2n\y (z - H + 2nl)* (z - H - 2nl)^ 

20 y^(x - X:) 2 ct^^(x - Xj) 2a/(x-Xj) 2o,^(x - Xj) 

e* '-i-e'' '+e '+e 


1\ 


f (HM 

I 1 OO 

' , + N' 


(II + 2nl)’ (ll-2nl)^ 


2a/(x.J 




/ 


(58) 


Since the true road length, S , will vary from the segmented length, 
we ensure conservation of pollution by scaling 


k 


2E Lj. 

i=l 


Qi = 


QS 


k 

Li 

i=l 
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Flirt liennorc, Qj may be varied along the route to account for high pollution 
points such as interchanges. 


ROADSIDE POLLUTION CONCENTRATION 


Using equation (58) as a specific example, it is apparent that the Qj’s are not 
yet specified. Other parameters, such as Xq, Xj, 0 ,^, Oy, etc. are available from other 
design estimates or tabulated data, but the Qj’s require further development. 

To begin the discussion, the emission rates in mass per unit length of travel from 
automobiles have been published in Reference 6. (These emission rates vary according 
to load.) Assume that under given conditions automobiles emit an average of Qj* grams 

per meter per vehicle. The only other dimensional quantity in equation (58) is the 
traffic density, N(t), which is measured in vehicles per second. Thus, equation (58) reads, 
dimensionally. 


[x] = [Q] [N] 


or 

Mass = jQi veh 
Length^ 'lime 


where: M is mass, L is length, and t is time. Assume a volume over which the 
emissions are stabilized. This volume is of height H, length L and expands in a 
direction perpendicular to the road at a rate u. The volume described is given by LjHu; 

where, as before, Lj is the length of the road which is considered. 

Thus, 


^ l^* _ 0 ^ 
LjHu Hu 


so that 


Mass -Time 
Length^ ■ Vehicles 


fQj] 


[Qj*] 

[H][u] 
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Obviously 


[xl 


Mass • Time Vehicles 
Length^ • Vehicles , . Time 


Mass 

Length* 


INPUT VARIABLES 


We shall now apply the theory which has been developed to sample data cases. 
The program used in the generation of this data is documented in Appendix A. 

Perhaps one of the most important sets of parameters given in the preceding are 
the diffusion coefficients Oy and a^- These quantities are fairly complicated functions 

of the distance from the source, x, and the meteorological category. From Reference 
7, we have the following tables. 


TABLE 1. DIFFUSION COEFFICIENT Oy, M 


Distance From 
Source, Meters 

Meteorological Categories 


A 

B 

C 

D 

E 

F 

10* 

22 

16 

12 

8 

6 

4 

10* 

210 

150 

105 

75 

52 

36 

10“ 

1700 

1300 

900 

600 

420 

360 

10* 

1 1000 

8500 

6300 

4100 

2800 

2000 


TABLE 2. DIFFUSION COEFFICIENT, a^, M 



A 

B 

c 

D 

E 

F 

10* 

14 

11 

7.6 

4.8 

. 3.6 

2.2 

10* 

500 

120 

70 

32 

24 

14 

10“ 

- 

- 

420 

140 

90 

46 

10* 

- 

- 

2100 

440 

170 

92 
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The meteorological categories can be characterized as follows [7] * . 


1 = A = extremely unstable conditions 

2 = B = moderately unstable conditions 

3 = C = slightly unstable conditions 

4 = D = neutral conditions (applicable to heavy overcast, night or day) 

5 = E = slightly stable conditions 

6 = F = moderately stable conditions. 

The categories A-F relate to surface wind speed, insolation, and percent of overcast 
as shown in the following table {?]. 


Surface Wind Speed 
Meters/Sec 

Daytime Insolation 

Thin Overcast or 
> '‘/s Cloudiness 

< V* Cloudiness 


Strong 

. 

Moderate 

Sliglit 



< 2 

A 

A-B 

B 



2 

A-B 

B 

C 

E 

F 

4 

B 

B-C 

C 

D 

E 

6 

C 

C-D 

D 

D 

D 

> 6 

C 

D 

D 

D 

D 


Some comments about the stability classes and corresponding Oy, values are 

in order. There exists a “G” category of stability which will not be considered here. 
This category requires a negative thermal radiation balance. The G category may thus 
occur only at niglit. The F category is also primarily a night-time reading, though F 
may extend beyond sunrise on the order of one hour. There is hesitancy on the part of 
some authors to use F (and G) class stability since an averaging process outside of the 
direct time period of applicability would yield false averages. References 8 and 9 treat 
F (and G) as common cases. 


1. For riuinerical convenience we shall relabel class A as 1 througli class F as 6. 
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SiiKv (lie period ol worst meteorological eonditioiis with respeet to pollution 
propagation does not oeeur at the same general time as does peak Irarfic loads, care 
must be taken that realistic combinations ace used. Even so. one of the following graphs 
will deal exclusively with F class stability in order to illustrate conditions which could 
exist under special circumstances. 

The numerical values of the a„, o_ parameters must also be used with care. 

\ ^ 

These are intended for fairly short range calculations* . Specifically, many of the values 
of Oy and given in Tables 1 and 2 are apparently the result of extrapolations rather 

than measured data. 

For this reason, we shall restrict the numerical use of the equations to a distance 
of one kilometer for x - Xq. 

Reference 10 contains valuable data concerning the calculation of the a values. 
In this source, the values of these parameters are shown as a functional equality involv- 
ing the wind equation angle, the vertical virtual distance, vertical diffusion, and the dis- 
tance over which rectilinear vertical expansion occurs downwind from an ideal point 
source. This approach could be used if one were attempting to isolate the exact set of 
meteorological parameters which correspond to worst case estimates. Even so. certain 
assumptions must be made with respect to the input variables to the equation unless 
massive quantities of statistics on meteorological patterns are available. For this reason 
it is convenient to use the values given in Tables 1 and 2 for values of x - x^ < 10* 
meters. 


There are other important parameters in the equatipns which have been developed 
besides the Oy, values. The parameter H, the height of stabilization of the center 

of the pollutant cloud, drastically affects the numerical results, H, or its equivalent, is 
almost inevitably included in all but the most elementary diffusion models. Specifically, 
rocket exhaust plumes during static firing have flow rates and heat contents which arc 
extremely large. This gives rise to very higli stabilization heights. Furthermore, stack 
gas emissions from power generation and industrial sources sometimes are designed to 
yield a maximum value of H within economic and materials constraints. In the case of 
automotive pollution, II is very much smaller than are the values for either rocket 
plumes or industrial stacks. 

An equation suitable for estimating H for automotive exhaust is given in Refer- 
ence 10. Assuming the radius of the exhaust pipe to be effectively zero we have, from 
110 ], 


2. Stephens, J. Brisco: Personal Communications. 



where: 


■ G is the buoyancy parameter, (3 gn/M^pCpT), 

? is the effective heat released, 

, p is the density of anbient air, 

Cp is the ratio of specific heats of air at constant pressure, 

r the temperature, 

, ^ the entrainment coefficient, 
g 90 

is defined at - — , and 
T" 9z 

g is the acceleration of gravity. 

^ n 

The quantity ^ is defined by the subsidiary relationship 

d6 _ dT ^ g 
, 9z 9z Cp 

j. Apparently fewer estimates have been given for the flovv rate and exhaust temper- 
rature of automobiles than have been given for certain other characteristics of these 
machines. Two sources which list the flow rate and temperatures of automotive engines 
are References 10 and 1 1. Reference 10 gives the following table, which is the more 
complete of the two sources. 

TABLE 3. ENGINE EXHAUST DATA 


i ; . j 

Driving Mode 

Idle 

Acceleration 

Cruise 

Deceleration 

Exhaust 

Flow 

C.F.M. 

5-25 

40^200 

25-60 

' 5-25 

m* /min 

0.14-5.66 

1.13-5.66 

0.71-1.70 

0.14-0.71 

Exhaust 
Gas Temp.* 

op 

150-300 

450-700 

400-600 

200-400 

°C 

65.6-149 

232-371 

204-316 

93-204 


* At the entrance to the muffler. 

























Additionally, a temperature profile for the atmospheric layer in the 1.2 to 7.1 
meter region is given in Reference 4. 

A very careful calculation would involve estimates of the percentage of time in 
each driving mode. But the calculations with limit cases indicate a rather high insensitivity 
to various inputs since H involves a fourth root. Specifically, although the temperatures 
listed above are almost surely much too high for the actual temperature of the gas emer- ' 
ging at the tailpipe, there seems no need for further precision. The estimates range from 
approximately one-half meter to about three and one-half meters, with two meters being 
a good number. This will be used in subsequent calculations. (The numerical range 
assumed for { for an H value of two meters is 0.11 to 0.16). 

In order to illustrate the primitive use of the diffusion relationship, a simple case 
will be presented, using a linear road which is orthogonal to the direction of the wind. 
Furthermore, a constant flow of 1000 vehicles per hour will be used in order that speci- 
fic traffic flow patterns do not muddle the overall picture. The length of the road will 
be specified at 2 kilometers. The mixing layer depth, 1, will be treated parametrically, 
and the road width given at 10 meters. 

Initial parameter variations will involve the distance from the road and the mete- 
orological categories. To obtain results of general application, x/Q will be reported 
rather than x and Q per se. The y and z coordinates used below are arbitrarily 
set as (0,2) respectively. 


NUMERICAL RESULTS 


The most obvious data which can be generated from the program concerns the 
reduction of pollution as one moves away from the road. This is plotted in Figure 7. 
Ihis data was generated under the assumption of class 6 stability - which is probably 
overly severe, but acts as a bounding case. Notice that for constant traffic density, the 
parameter u enters hyperbolically, so that only one value of ii needs to be actually 
programmed. 

For a given wind speed, it is also of interest to illustrate the variation of x/Q 
with distance for varying stability classes. This is shown in Figure 8. This Figure also 
assumes a mixing layer depth of 1000 meters. The figure shows the importance of 
stability class variation with respect to pollution concentration. For a very unstable 
atmosphere (stability class 1 or A), pollutants are quickly dispersed due ito upward 
ventilation. As the atmosphere becomes progressively more stable, less and less vertical 
mixing occurs, resulting in higher concentrations near the ground. ! 
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1.0 



Figure 7. Variation of x/Q versus distance from line source, with wind 
speed as a parameter (stability class 6(F)-mixing layer depth 1000m). 

Rather than dealing directly with varying x/Q. it is also of interest to fix the 
value of x/Q via variation of wind speed. That is, for a given jdistance there is some 
wind speed which corresponds to a value of x/Q ” • Since distance and wind speed, 

at which \/Q “ , depends upon the stability class and the mixing layer depth, these 

additional variables must also be reckoned with. In Figure 9, the mixing layer depth is 
1000 meters and the stability categories are shown parametrically. Thus, for a given 
stability class and a mixing layer depth of 1000 meters, we can find the distance at which 
x/Q is equal to e ^ for a given wind speed. 

In Figure 10 we have the same plot as Figure 9, excepting the fact that the 
mixing layer depth is only 100 meters. It is interesting to note that even this drastic 
decreases modifies the results only for stability class 1 and then only for distances be- 
yond several hundred meters. The effect would eventually modify all stability classes 
at further distances, of course. 
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Figure 8. Variation of x/Qversus distance from line source, with 
stability class as a parameter (wind speed 1 m/sec — 
mixing layer depth 1000 meters). 

Figure 1 1 indicates the result of a further reduction in mixing layer depth. The 
case illustrated has a mixing layer depth of only 10 meters - a low but physically realistic 
value. The result is obvious for all stability classes - namely, the pollutants are trapped 
near the ground and concentration becomes almost independent of distance. 

The next three graphs. Figures 12, 13, and 14 are cross plots of the data shown 
in the preceeding three graphs. In each case, the wind speed necessary to achieve the 
value of x/Q* = at the parametric distances is plotted against mixing layer depth. 
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i ‘ . 

Figure 9. Distance from line source at which x/Q = 1/e^ versus 
wind speed, with stability class as a parameter (mixing layer depth ,1000 meters). , 
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DISTANCE FROM LINE SOURCE, METERS 
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Figure 10. Distance from line source at which x/Q = e ^ versus wind speed, 
with stability class as a parameter (mixing layer depth 100 meters). 
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WIND SPEED, METERS/SECOND 


Figure 1 1. Distance from line source at which x/Q = e ^ versus wind speed, 
with stability class as a parameter (mixing layer depth 10 meters). 
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Figure 13. Wind speed required to achieve 
stability class (mixing 


34 


5 


; parametric distance V( 
leters). 









1U 

9 


8 

7 

6 




630 1000 


(COINCIDENT) 


NOTE: 


INTERNAL NUI 
(METERS) TO 


MBERS INOICATp 
LINE SOURCE 


DISTANCE 


1 2 3 4 5 6 

STABILITY CLASS 


Figure 14. Wind speed necessary to achieve x/Q = e'* at the parametric 
distance versus stability class (mixing layer depth 10 meters). 

The three figures differ in having mixing depths of 1000, 100, and 10 meters, respectively. 
The representation of the data in this form is particularly useful if a statistical study of 
probable weather pattern is to be used to predict potential pollution levels. 

In the above mentioned charts, the mixing layer depths were never taken above 
1000 meters. The reason for this limitation is that 1000 meters is an asymptotic value 
insofar as diffusion over the x coordinate range of 1000 meters is concerned. In other 
words, values of mixing depth above 1 000 meters give no significant variation from a 
mixing layer depth of 1000 meters. 

This point is further illustrated in Figure 15. The wind speed necessary to achieve 
a x/Q of e’* at 1000 meters versus mixing layer depth is illustrated in this case, with 
stability class as a parameter. From this graph it is apparent that lower mixing ceilings 
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WIND SI»EED. METERS/SECOND 


!• 



10 100 1,000 10,000 


MIXING LAYER DEPTH, METERS 

Figure 15. Wind speed necessary to achieve x/Q = e ^ at 1000 
meters versus mixing layer. 

most strongly affects the most unstable atmospheric conditions. This may be intuitively 
understood as follows. The basic equations for describing the mixing lid reflections is 
the method of infinite imaging, equation (16) and the modifications that follow from it. 
Under extremely unstable conditions (say, 1 or A conditions) then a cell of pollutant 
would tend to rise very rapidly, and thus be “reflected" at an earlier point than would 
a cell of pollution under more stable conditions. Thus, Figure 15 corresponds to physical 
intuition. 
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CONCLUSIONS 


It has been shown that the standard diffusion equation may be separated to 
yield solutions which interface with known geometric and meteorological constraints. 
Specific weather conditions were assumed to illustrate the application of the theory. 

It was seen that under the assumption of constant traffic density on a linear 
road, the worst weather conditions involve a rather stable atmosphere and low wind 
speed. The depth of the mixing layer seems relatively unimportant of the first kilometer 
from the line source, unless very low mixing depths are encountered. 

Under the very special conditions of small mixing depths, the most stable 
atmosphere may, at adequate distance, actually produce a lower concentration of pollu- 
tants than an unstable atmosphere. 
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APPENDIX A 


BASIC COMPUTER PROGRAM 
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SUBROUTINE LG I NT ( L AGRANGI AN INTERPOLATION) JRA 12/12/72 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

C SUBROUTINE LGINT PERFORMS LAGRANGIAN INTERPOLATION USING A C 
C SPECIFIED NUMBER OF ' THE DATA POINTS IN THE GIVEN TABLE. C 
C C 
C ' C 
C SUBROUTINE'S ARGUMENTS ARE... C 

C c 

C XTABL — ARRAY OF INDEPENDENT VALUES OF- DATA TABLE. C 
C C 
C YTABL — ARRAY OF DEPENDENT VALUES OF DATA TABLE* VALUES C 
C ARE IN DIRECT CORRESPONDENCE WITH XTABL VALUES. C 
C , C 
C NTABL — LENGTH OF XTABL AND YTABL ARRAYS. C 
C C 
C X —INDEPENDENT VALUE FOR WHICH A CORRESPONDING C 
C DEPENDENT VALUE IS TO BE INTERPOLATED. C 
C C 
C Y — THE interpolated DEPENDENT VALUE. C 
C C 
C NPTS — NUMBER OF POINTS FROM THE TABLE TO BE USED IN THE C 
c interpolation. C 
c C 
c ierr — error indicator ierr -1 interpolation performed c 
c correctly c 
c IERR*2 interpolation INHIBITED C 
C BECAUSE X OUTSIDE TABLE* C 
C NPTS LESS THAN 1* OR NPTS C 
C GREATER THAN NTABL. C 
C C 
C SUBROUTINE'S IMPORTANT VARIABLES... C 
C C 
C NONE c 
C C 
C LIMITATIONS... C 
C C 
C values of XTABL ARRAY MUST BE STRICTLY INCREASING. C 
C c 

C X GREATER THAN OR EQUAL TO XTABL! 1) AND LESS THAN OR EQUAL C 
C TO XTABL<NTABL). C 
C C 
C NPTS GREATER THAN OR EQUAL TO 1 AND LESS THAN OR EQUAL C 
C TO NTABL. C 
C C 
C NOTES. •• C 
C C 
C LGINT IS LIMITED TO INTERPOLATION* IT WILL NOT . EXTRAPOLATE. C 
C C 
C ONE POINT INTERPOLATION RETURNS THE YTABL VALUE WITH C 
C THE CORRESPONDING XTABL VALUE THAT IS CLOSEST TO X. C 
C C 
C TWO POINT INTERPOLATION IS EQUIVALENT TO LINEAR C 
C INTERPOLATION USING THE TWO POINTS. C 
C C 
C THE TABLE VALUES USED IN THE INTERPOLATION ARE CENTERED C 
C AROUND THE POINT OF INTERPOLATION WHEN POSSIBLE* HOWEVER* C 
C WHEN THE INTERPOLATION POINT IS CLOSE TO EITHER END OF C 
C THE TABLE* THE TABLE VALUES USED ARE THAT END OF C 
C THE TABLE. C 
C C 
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c C PROGRAMMER— JAMES R* ALEXANDER C 

c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE LGI NT { XTABU • YTARL *NT ABL »X t Y *NPTS * I ERR ) 

DIMENSION XTABL( 1 ) »YTABL( 1 ) 
lERR-X 
C 

C ERROR CHEQKS 

C 

IF(NPTS)20t20«10 
10 IF(NPTS-NTABU30*30*20 
20 IERR-2 
RETURN 
C 

C INTERPOLATION INHIBITED— VALUE OF NPTS OUT OF PERMISSIBLE RANGE 

C LOCATING PORTION OF TABLE TO BE USED IN INTERPOLATION 

C 

30 1-1 

I2»2*(NTABL-1) 

NC«NTABL 

11- l 

40 IF(X-XTABL( I ))60*60»70 
50 Y«YTABLm 
RETURN 
60 I2>I 
NC«I 

GO TO 80 
70 Il-I 

80 I*( I2-Il+l)/2+Il 
I2>NC 

IF(I2-I1-1)90.100»h0 
90 IERR-2 
RETURN 
C 

C interpolation INHIBITED— value OF X OUTSIDE RANGE OF XTABL VALUES 
C 

100 NC-I2 

IF(XTABL( I2)-fXTA6L( I1)*2.0*X)120»120*110 
110 NC-Il 

120 1 1« I l-NPTS/2fl- ( NPTS-NPTS/2*2 > » ( 1 1/NC) 

12- I1+NPTS-1 

IF( II )130*130»140 
130 Il-l 

I2-NPTS 
GO TO 160 

140 IF{12-NTABLI160»160*150 
150 Il-NTABL-NPTS-fl 
I2-NTABL 
C 

C PERFORM INTERPOLATION 
C 

160 Y-0.0 

DO 190 I-11.I2 
YY-1.0 


DO 180 J-I1«I2 
IF(I-J)170.180.170 

170 YY-YY*(X-XTABL( J» )/(XTA8L< II -XTABL tJ I ) 

180 CONTINUE 

190 Y«Y+YY*YTABL< I ) 

RETURN 

END 
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FUNCTION SIGY(X.VIS) 


FUNCTION SIGY DETERMINES THE HORIZONTAL DISPERSION COEFFICIENT 
TABLE OF values TAKEN FROMi WORKBOOK OF ATMOSPHERIC DISPERSION 

ESTIMATES* 0. BRUCE TURNER* PUBLIC 
HEALTH SERVICE* 1970— 'AVAILABLE AS 
PB-191A62 OR N71-31626* 

SEE FIGURE 3-2> 

X —distance DOWNWIND 

VIS —STABILITY CLASS 1>A**..»6>F* A NON INTEGER VALUE' RESULTS 
IN linear INTERPOLATION IN THE LOGX-LOGSIGY SPACE BETWEEN 
THE TWO NEAREST CLASSES 

DIMENSION FY(16*6) 

DIMENSION FI (16) *F2 ( 16 ) *F3 ( 16 ) *F4 1 16 ) *F9 ( 16 > *F6 ( 16 1 
EOL I VALENCE (FY(l)»Fl(in*(FY(17)*F2 (!))*( FY( 33 )*F3(1)) 
EQUIVALENCE ( FY ( 49 ) *F4 ( 1 ) ) * ( FY(65 I *FS 1 1 1) * ( FYIBl ) *F6 1 1 1) 

DATA F1/27.6*41.6*62.0*95*0*140«0*215.0*32S*0*470*0*715.0* 

11040. 0*1590.0 *2300* 0*3350.0 *4900*0 *7250.0 *11000*0/ 

DATA F2/19.5*29.0*44.5*6e*0*100*0«15e.0»236* 0*360*0 *550*0 *800*0* 
11190. 0*1800. 0*2600*0 *3850. 0*5580*0 *8200*0/ 

DATA F3/12. 5 *19. 2 *34.0 *45. 5 *68*0 *105.0 *162*0 *242*0*370*0* 

1550. 0*840. 0*1280. 0*1880. 0*2800. 0*4150*0 *6100*0/ 

DATA F4/8. 1*12. 3 <19. 0*26. 5 *44. 0*68*0* 108 *0*160*0*247*0*370.0* 
1560.0*840*0*1220.0*1820*0*2700*0*4100*0/ 

DATA F5/6. 1*9. 3 *14.2 *21 *8 *33.0 *5 1*0 *78*0* 118*0*180*0*275*0* 
1410.0*630*0*830*0*1380*0*2050*0*3050*0/ 

DATA F6/4. 1*6.2 *9.3 *14*8 *22*0 *34*0 *53*6*79*0 *122*0 *180*0* 
1275.0*420.0*620.0*940.0*1380*0*2050*0/ 

A(X)«0*43429»AL06(X) 

N = 1 

IF(X-1*0)10*20.20 
10 SIGY*0*0 
RETURN 

20 IF(X-100000. 0)30*30*10 
30 I1*VIS+0*1 
I2-11+1-I1/6 
IF(I1 )10*10*40 
40 IF( 11*6)50*50*10 
50 IF(X-160*0)70*70*60 
60 N>l*0001-t>(A(X)*2*0)«5*0 . 

70 Yl«A(FY(N*il) )+A(FY(N*12)/FY(N*Il) )*( VIS-FLOAT ( I 1 ) ) 

Y2»A(FY(N+1*I1 ) )+A(FY(N+l*I2)/FY(N4l*Hn*(VIS-FL0AT<Il) ) 
SIGY«10.0**( Yl+I Y2-Y1 ) /0*2«( A(X )-2*0-(N-l )»0*2) ) 

RETURN 

END 
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onnnnnnnnnno 


FUNCTION SIGZ(X*V1S) 


FUNCTION SI6Z DETERMINES THE VERTICAL DISPERSION COEFFICIENT 
TABLE OF VALUES TAKEN FROM* WORKBOOK OF ATMOSPHERIC DISPERSION 

ESTIMATES. 0. BRUCE TURNER* PUBLIC 
HEALTH SERVICE* 1970— AVAILABLE AS 
PB-191<f82 OR N71-31626* 

SEE FIGURE 3-3. 

X — DISTANCE DOWNWIND 

VIS —STABILITY CLASS 1*A*...*6»F. A NON INTEGER VALUE RESULTS 
IN LINEAR INTERPOLATION IN THE LOGX-LOGSIGZ SPACE BETWEEN 
THE TWO NEAREST CLASSES 

dimension FZ(16*6) 

DIMENSION FI (16) *F2 ( 16 ) *F3 ( 16 ) *FA( 16 ) *F5 ( 16 ) *F6 ( 16 ) 

EQUIVALENCE (FZ ( 1 ) »F1 ( 1 ) ) * < FZ ( 17) *F2< 1 ) ) . (FZ ( 33 ) *F3 ( 1 ) ) 
EQUIVALENCE ( FZ( 49 ) *F4 ( 1 ) ) * ( FZ ( 65 ) *FS( 1 ) ) * I FZIBl ) *F6 ( !•) ) 

DATA FI /14.0 *23. 0*37 *0*72*0 *170.0 *465.0* 1220 *0*3100.0 *8600*0* 
121500. 0*60000.0* 140000 .0.31500Q.0* 740900*0 *1650000. 0*4200000*0 
DATA F2/11 .0*16*5 *25.5 *41*0 *65.0* 110*0 *185.0 *300.0*500*0* 
1820.0 *1380.0 *2250. 0*3750*0 *6200.0* 10200*0 *17000*0/ 

DATA F3/7.4* 11 *5 *17.0 *26*3 *39.5 *61 *0*95*0 *142 *0*220.0 *330.0* 
1510.0*780.0*1150*0*1800*0*2650*0*3900*0/ 

DATA F4/4. 7 *7* 0*10* 3 *15. 8*22 *2 *32 .0*43*8 *58*0 *77*0 *101*0* 
1138.0*180.0*230.0*295*0*365.0*460*0/ 

DATA F5/3. 6 *5. 3 *7 *7 *11 *0*15.3 *21.2 *29*0*36*5 *50 .0*62*0*79*0* 
1100.0*120.0*143*0*163*0*183*0/ 

DATA F6/2 *3 *3*4 *4 *9 *7. 1*10*2 *14*0* 19*0 *24*0 *31 *0*37*5 *46*5* 
156.0*64.0*73.0*84*0*94*0/ 

A(X)»0.43429*ALOG(X) 

N»1 

IF(X-1. 0)10*20*20 
10 SIGZ>0*0 
RETURN 

20 IF(X-100000.0)30*30*10 
30 I1>V1S40*1 
I2*>Il4-l-Il/6 
1F(I1)10*10»40 
40 IF(I1-6)50*50*10 
50 IF(X-160*0)70*70*60 
60 N»l*0001+(A(X)-2*0)*5*0 

70 Y1«A(FZ(N*I1) )+A(F2(N*I2)/FZ(N*Il) )*( VIS-FLOAT ( I 1 ) ) 

Y2«A(FZ(N+1*I1 ) )+A(FZtN+l*I2)/FZ(N+l*Il) )*(VIS-FLOATI ID ) 
SIGZ-10«0**(Yl+( Y2-Y1 )/0«2*(A(X)-2*0-(N-l)*0*2) ) 

RETURN 

END 
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FUNCTION CI:IIF3«FX»Z»HIL»M) 

C 

C FUNCTION CHIF3 EVALUATES THE 2-X DEPENDENCY OF CHI 
C FX —THE VERTICAL DISPERSION COEFFICIENT AT X 

C Z — Z COORDINATE 

C HIL —HEIGHT OF INVERSION LAYER 

C H —SOURCE HEIGHT 

C SERIES IS SUMMED UNTIL NTH TERM CONTRIBUTES LESS THAN 0*01 i>ERCENT 
C TO the N terms of THE SERIS. 

C 

C»0.5/<FX*FX> 

R«Zi-H 

S-Z-H 

F«EXPI-C*R»R)+EXP(-C*S*S» 

N*1 

10 0«FL0AT(2*N»*HIL 

FN*EXP<-C*(S-0»*«S-0) »+EXPI-C*(R+Q)*IR+0» UEXP(-C*(R-0)»(R-0) I 
l+£XP<-C*(S+0»*(S+0) » 

F-F+FN 

N«N+l 

IF(FN/F-0*0001)20*10*10 
20 CHIF3*F 
RETURN 
END 


FUNCTION CHIFA(HX»Y»YO) 

C 

C FUNCTION CHIF4 EVALUATES THE Y-X DEPENDENCY OF CHI 
C HX —THE horizontal DISPERSION COEFFICIENT AT X 
C Y — Y COORDINATE 

C YO —ROAD SEGEMENT HALF LENGTH 

C 
C 

C REAL ERROR FUNCTION APPROXIMATION 

C ERF(X) « ERFA(X) E* X GREATER THAN OR EQUAL TO ZERO 
C ABS(E) LESS THAN OR EQUAL TO 1.5E-07 

C REFERENCE* ABRAMOWITZ* M* AND STEGUN* I* A* HANDBOOK OF 
C MATHEMATICAL FUNCTIONSt NATIONAL BUREAU OF STANDARDS* 

C APPLIED MATHEMATICS SERI£S*SS* JUNE 1964* PAGE 299* 

C EQUATION 7.1.26* 

C 

S(X)»I.0/( 1.0+0 *3275911*X) 

ERFA(X)»l.0-(S(X)»(0*254e29592+S(X)*«-0*284496736+S«X)* 

1 ( l*421413741+S(X)«(-l*4S31S2027+S(X)*l*06140»429n Ml* 

2 EXP(-X*X» 

C 

C 

C-HX»1. 41421356 

Yl-YO-Y 

Y2-Y0+Y 

CH I F4»S1 GN ( ERF A I ABS < Y 1 1 /C ) * Y1 1 +S IGN I ERFA < ABS I Y2 » /C ) • Y2 1 

RETURN 

END 



** PROGRAM CHITR(PRE0ICT10N OF STRAIGHT ROAD AIR POLLUTION! 07/28/73 
DIMENSION TT(50) .PHI <50) 

IRO«2 
I PR-3 

C SIGZM IS The vertical dispersion assumed TO HAVE OCCUREO 
C FROM VEHICLE INDUCED MIXING* 

SI62M-2.4 

C READ TRAFFIC TIME DENSITY! TIME IN UNITS USED BELOW* VEHICLES/SEC I 
C DENSITY FUNCTION IS USED CYCLICALLY 

READCIRDtl) N*(TT(I)*PHI(1)*I-1*N) 

1 FORMAT! I5/I2F0. 2) > 

TTT-TT!N)-TT!1» 

C READ EMMISSION/METER-VEH* ROAD HALF LENGTH*S0URCE HEIGHT* INVERSION 
C LAYER HEIGHT. ROAD HALF WIDTH* WIND VELOCITY* AND STABILITY CLASS 
C ! 1 «0— A* • • * *6*0— F I 

READ!IRD*2) Q*YO*H*HIL*XM*U*VIS 

2 FORMAT! 10F8. 2 » 

C CALCULATE X BIAS FOR INITIAL DISPERSION 
Y1>SIGZ!1.0*VIS) 

Y2-SIGZ!100.0«VIS) 

XB«10*0*«!2.0»AL06!SIGZM/Yl)/ALOG!Y2/Yin 

SIGYO>SIGY!XM-«-XB*VIS) 

SIGZ0>SIGZ!XM-4-XB*VISt 

C1-1.0/CHIF3!SIGZ0*0.0*HIL*H) 

C2«1.0/CHIF4!SIG,YO.O.O*YO) 

C-0/H/U*SIGZ0*Cl*C2 
C READ NUMBER OF PREDICTIONS 
READ!IRO*7) NPRED 

7 FORMAT! 151 

DO 100 IPRED-1*NPRED 

C READ SPATIAL AND TEMPORAL COORDINATES OF PREDICTION POINT 
REAO!IRD*2i X*Y*Z*T 
IF!X-XMJ75*40«40 

40 IF!Z»75*50*50 

50 IF!Z-HIL)80*80*75 

75 WRITE11PR*5I 

5 FORMAT! 1H1/10X*26HINPUT COORDINATES IN ERROR) 

GO TO 100 

80 SIGYX-SIGY!X-*XB*V1S> 

SIGZX-SIGZ!X*-XB*V1S) 

Fl-l.O/SIGZX 
TP-T-!X-XM) /U 
M«ABS!TP/TTT) 

1F!TP-TT!1) 182*83*83 

82 TP«TP+TTT*!M+1) 

83 IF!TP-TT!N) )85. 85*84 

84 TP-TP-TTT*M 

85 CALL LGINT!TT*PH1 .N*TP.F2.2»IERR) 

GO T0!95*90) *1ERR 

90 WRITE!lPR*6) TP 

6 FORMAT!1H1/10X*19HINTERPOLATION ERROR/10X*F8.2 ) 

GO TO 100 

95 F3-CHIF3!S1GZX*Z*H1L*H) 

F4-CH1F4IS1GYX*Y*Y0) 

CH1-C*F1*F2*F3*F4 
WRrTCIlPR*3) X*Y*Z*T*CHI 

3 F0RMAT!lHl/10X*4HCHt!*F8*2*lH**FS*2*lH**F8*2*lH**F8*2*4H) > * 
1E12.5) . 

WR.TE!IPR*4) O*Y0*H*HIL*XM*U*VIS*XB*SIGY0*StGZ0*SlGYX*SlGZX* 
1C1*C2.C*F1»F2*F3.F4 

4 FORMAT! //10X.4HQ - *E 12 .5/lOX . 5HY0 - .E12.5/ 
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U0X»4HH - tE12tS/10X»6HHlL > E12*»/10X*$HXM ■ tE12«»/10X(4HU • t 
2E12*&/10X»6HV1S > *E12*5/10X*5MX8 • tE12*9/10X »aHSIGYO > *E12*»/ 
310X*8HSIGZ0 - *E12*3/10X*8HSIGYX « tE12«5/10X*8HSlGZX > t€12*5/ 
410XfSHCl > *E12.5/10X»5HC2 > »E12*»/10Xt«HC • f E12*9/10Xr8HFl > • 
5E12.5/10X»8HF2 > «E12*5/10X(8HF3 > tE12.S/10X»5HF4 > »E12«S) 

100 CONTINUE 
STOP 
END 



APPENDIX B. EVALUATION OF CONTOUR INTERVAL 


As previously mentioned, it is usually necessary to express pollution concentrations 
from a road when the wind vector is not orthogonal to the road. The specific integral for 
this task is shown on page 20. 

In order to save computer time in the implementation of this case, it is often 
necessary to approximate the given integral analytically. To this end we proceed as 
follows. 


If we agree to evaluate the integral for very large fi and at z = H we have 


Vo tan 

X (X, y. H, t) = ^ (Xo) Oy (Xq) J 

-Vo tan < 1 / 


. (y - X* tan_i^/_ 
g 2 0y^(x-x*)‘ y b 




- 

Oy (X - X*) Oz (x - X*) 




dx* (B-1) 


We now assume that 




which will be valid for short segments. Then (B-1) becomes 


(y - X* tan i//y 


yo tan i// 2 Oy^ (x - x*) 

X (X. y, H, t) = >- (Xo) Oy (Xo) N (t) J — rr dx" 


-yo tan 


Oy (x - x*) O2 (x - X*) 


For convenience we now write. 


1 = 


I 


(y - X* tan 1//)^ 

2 Oy^ (x - X*) 
e ^ 

Oy (x - X*) O2 (x - X*) 


dx* 


(B-2) 
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Then, 


y - X* tan 1// = y - X tan + X tan (// - X* tan = (y - X tan i//) + (x - 

X*) tan \l/ 

Abbreviating, 


(x - X*) = r 

(B-3) 

y - X tan \f/ = r} , 

(B-4) 

we have 


2 

(77 + r tan \l / ) 

f 2 Oy' (r) 

I-- « dr . 

J Oy (r) (r) 


We now approximate the reciprocal Oy in the exponential as 


joy (r)J'‘ = (6 + e r) 

(B-5) 


and operate on the resulting exponent as follows: 

(t? + r tan i//)^ (6 + e r)^ = Fs rj + r (5 tan i// + e 17) + (e tan \J/) r^ 


Substitute 


2 e tan 


(B-6) 


so that 


(Tj + r tan \jj)^ (6 + e r)^ = (e tan \p) s^ 


t:ni ^ ^ 
4 € tan ^ 
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Then 


(5 tan i// - e t?) 
j ^ g 4 e tan 


f Q-{etmyp)s^ 

J a^{s + ix) a, (s 


(s *:• m) 


ds * , 


with 


_ (6 tan \l / -e T}) 
^ 4 e tan 


Expand the Oy and values that occur in the demoninator as 


ay(S + p) j-Q 


= Yj V * (s + p)* 


= y* V:** (s + ju)’ 

a^(s + M) > 


so that 


-} = y V. (s + ju)* 

Oy (S + P) Oy (S+ fi) 


(B-7) 


(B-8) 


(B-9) 


where the ^j’s are defined as obvious combinations of i/j* and . Then 


(5 tan ^ - € rj)^ 
I _ g 4 e tan t// 


Z / 


-(etan «//)s" (s + p)'ds 


Finally, setting 


y/ e tan \{/' s = 


S = CO 


(B-10) 
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so that 


(6 tan \j/ - e 17)^ 


I = - e 


4 e tan 


S “ — ~ — i+1 I ® ^ >/e tan ^ 

i=0 J 


(e tan 


e tan 1^)* dc^ . 

(B- 11 ) 


The integral may now be evaluated from expressions contained in standard 
.references such as Reference 12 . From that source, 


J e'*^ d X = erf x 
j 2 


J 


xe'^^d * " ■ ^ 


r x” e*’^^ d X = - — x*'"^ e"*^ +-^ / d x 
J 2 2 •' 


The coefficients which occur in the expansion of the binominal are denoted as 
Do = (e tan (»>o + »^iM + 

Di = (e tan (t', + Iv^ix St'a/u* + 4*'4/u’ + + 6Vf,ti^) 

Dj = (e tan {y^ + ‘iv-^n + 6 v^n^ + lOi^sP® + 

D3 = (e tan (1^3 + 4P4H + ICh'jp* + 20vo/i’) 

D4 = (e tan \(/y^/^ (v^ + SvjM + ISmsM*) 

Ds = (e tan ^y^ (Vs + bi'oM) . and 

Do = (e tan r'o (B-12) 
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or by the convenient binomial coefficient matrix as 


Do 

(e 

tan 

^)vr 


1, 

M> 



M^, M®, 

M® 


^0 

D, 

(e 

tan 



0, 

1, 

2/i, 

3m* 

, 4m* , 5m^ , 

6m® 



D, 

(6 

tan 



0. 

0, 

1, 

3m, 

6m*, 10m* 

, 15m“ 


^2 

Da 

(e 

tan 


= 

0, 

0, 

0, 

1, 

4m, 10m* 

,20m* 



D4 

(e 

tan 



0, 

0, 

0, 

0, 

1, 

5m, 

15m* 



Ds 

(e 

tan 

^)3 


0, 

0, 

0, 

0, 

0, 

1, 

6m 



D* 

(e 

tan 



0, 

0, 

0, 

0, 

0, 

0, 

1 


^6 


Then, 

(5 tan \p-€r))^ 

I = . e ^ ^ 

1 2 3 1 S 

- 2 [(Di + Da + SDj) + (D^ + - D4 + - D^) co + (D3 + 2 Ds)o>^ 

+ (D 4 + I D6)o) 3 + DsCo“ + Dsco*] I (B-13) 


16 


(8 Do + 4D2 + 6D4 + 15Do) erf co 


The choice of coordinate system (at the center of the given segment) guarantees 
symmetric integral limits of the variable x*. But it should be remembered that a linear 
translation of this variable was made, so it cannot be argued, in general, that the even 
terms in the last expression vanish. If we reintroduce the initial variables we have 




(x - X*) - 


[ 


8 tan 4 / + e(y - x tan 
2 € tan \jj 


]| 


\/e tan 


Thus, 
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‘^U*=yu lan \i/ 


X - yo tan 


[b tan tp + e(y - x tan i//) |1 , ■ 

2eun* - JJ V5T555T- 


(B-14) 


and 


^lx*=-yo tan ./. = ^ tan P 


[ 


6 tan + e(y - x tan 

(B-15) 


2 e tan p 


1 


Finally abbreviating, 


Eo=^ (8D„ + 4D2 + 6D4 + ISD^) 
1 6 


E, =-- (D, + D3 + 3 Ds) 


Ej = ^ (D2 + ^ D4 T- — L/6 


3 15 ^ , 

- D4 + ^ De) 


E3 = 2 (D3 + 2Ds) 
E4 =-!- (D4 +1 De) 


2 ' ’ 2 
1 


Es = 2 Ds 


, and 


E« = - D. 


We have 


I = + e 


(5 tan \p - €T})^ I 

4 e tan \p | q-oj^ [E, + Ejco + £ 3 ^^ + E^co^ + EjW^ + E(,a)*l 


- Eq erf CO > I 


(O, 

CO2 


(B-16) 


so that 


52 


I I I I I Mil 


I ■ II I I 



(6 tan \p - er))^ . 

X = 5 (X.) „y (x„) e‘ 4V tanir~ I 

+ E2 (coie'^^i^ -co2e'^2^) + E3 (o),^e‘^i^ -co2^e"^2^) + £4 

- 0)2^ e‘^2^) + Es (co,‘’e'^i^ -a)2^e'^2^) + Eg (a>i®e'^i* -W2^e'^2^) 

- Eo (erf CO, -erf oj 2 )l N(t) (B-17) 


It would be very convenient if equation (B-17) could be inverted to yield y as 
a function of x< since this would allow noninterpolative plotting of constant x isopletiis. 
Such an inversion is not a simple matter, however. 
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